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Abstract 

Fractional Brownian motion (fBm) structures are used to investigate the dependency of column density 
variance (uf DN ) in the turbulent interstellar medium on the variance of three-dimensional density (crj^ ) and 
the power-law slope of the density power spectrum. We provide quantitative expressions to infer the three- 
dimensional density variance, which is not directly observable, from the observable column density variance 
and spectral slope. We also investigate the relationship between the column density variance and sonic Mach 
number (M s ) in the hydrodynamic (HD) regime by assuming the spectral slope and density variance as functions 
of sonic Mach number, as obtained from the HD turbulence simulations. They are related by the expression 
°fn a< = ^°fn p = ^ m ( 1 + b 2 M^), suggested by Burkhart & Lazarian for the magneto-hydrodynamic (MHD) case. 
The proportional constant A varies from « 0.2 to ss 0.4 in the HD regime as the turbulence forcing parameter 
b increases from 1/3 (purely solenoidal forcing) to 1 (purely compressive forcing). It is also discussed that the 
parameter A is lowered in the presence of a magnetic field. 

Subject headings: ISM: structure — ISM: clouds — turbulence 



1. INTRODUCTION 

A variety of observations and simulations have shown 
that the density structures of the interstellar medium (ISM) 
are scale-free, hierarchical, and fractal (e.g., Stutzki et al. 
1998; Elmegreen & Scalo 2004; Burkhart et al. 2012). The 
supersonic, compressible turbulence is likely responsible for 
the complex and hierarchical density structures observed in 
the ISM. The density fluctuation due to the supersonic tur- 
bulence plays a crucial role in models of star formation 
rates and initial mass functions (Mac Low & Klessen 2004; 
Krumholz & McKee 2005; Hennebelle & Chabrier 2008; 
Federrath & Klessen 2012). Hence, the statistical proper- 
ties of density structures in the turbulent ISM have been 
extensively studied in hydrodynamic (HD) and magneto- 
hydrodynamic (MHD) simulations. 

It is now well known that the probability distribution 
functions (PDFs) of the three-dimensional (3D) densities 
and column densities of the turbulent ISM are close to 
lognormal (Vazquez-Semadeni 1994; Nordlund & Padoan 
1999; Klessen 2000; Ostriker et al. 2001; Wada & Norman 
2001; Burkhart & Lazarian 2012). The standard deviation 
of density (c p/ / po ) increases with sonic Mach number, M s 
(Nordlund & Padoan 1999; Ostriker et al. 2001). This trend 
is expressed by cr p / po = bM s , where p/po is the density 
{p) normalized by the mean density (po) and b is a con- 
stant of proportionality, known as the turbulence forcing 
parameter. Values of b ranging from « 0.3 to w 1.0 
have been suggested in numerical simulations (Padoan et al. 
1997; Passot & Vazquez-Semadeni 1998; Kritsuk et al. 2007; 
Beetz et al. 2008). Federrath et al. (2008, 2009, 2010) have 
found that for the same M s , compressive (dilatational) forcing 
leads to much larger density variance compared to solenoidal 
(rotational) forcing. They showed that b depends on the type 
of the turbulence forcing varying from k 1/3 for solenoidal 
forcing to « 1 for compressive forcing. 
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The density power spectrum is also a useful tool to char- 
acterize the turbulent ISM. The power spectra of turbu- 
lent clouds are power laws in form (Crovisier & Dickey 
1983; Despande et al. 2000; Stanimirovic & Lazarian 2001; 
Padoan et al. 2004). Kim & Ryu (2005) have investigated the 
dependency of the power-law slope of the density power spec- 
trum on M s and found that the density spectrum becomes 
gradually shallower as M s increases in HD turbulent me- 
dia. Kowal et al. (2007) and Burkhart et al. (2010) confirmed 
the flattening of the density spectra with M s in MHD cases. 
Federrath et al. (2009) showed that compressive forcing leads 
to significantly steeper density spectra than solenoidal forc- 
ing. 

We note that the spectral slope of the 3D density power 
spectrum is the same as that of the projected column density 
power spectrum (e.g., Stutzki et al. 1998; Padoan et al. 2004) 
and thus can be directly extracted from observations. On the 
other hand, the 3D density is not an observationally accessi- 
ble quantity, but instead the column density, i.e., the integral 
of the density along a line of sight, is observable. There is 
a great deal of observational data, which traces the ISM col- 
umn density, from surveys. Some examples are the GALFA 
(Peek et al. 201 1) and GASKAP (Dickey et al. 2012) H I sur- 
veys and Columbia-Cfa (Dame et al. 2001) and ThrUMMS 3 
molecular gas surveys. The best way to infer the variance of 
the 3D density field from observations is obviously to com- 
pare the observational data with a large number of turbulence 
simulations. However, instead of performing expensive simu- 
lations, a simpler approach to derive the 3D density variance 
would be of great interest. Padoan et al. (1997) generated ran- 
dom density distributions with various density variances and 
spectral slopes and compared the resulting synthesized ex- 
tinction maps with the observed extinction data. Brunt et al. 
(2010) developed an interesting technique for calculating the 
variance of a 3D density field using only information con- 
tained in the two dimensional projection. The technique has 
been applied to the Taurus molecular cloud (Brunt 2010). 
More recently, Burkhart & Lazarian (2012) found a column 

3 http://www.astro.ufl.edu/~peterb/research/thrumms/ 
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Figure 1. Examples density power spectra with the 128 3 , 256 3 , and 512 3 
resolutions, which were obtained for (7 g , a\ np ) = (3.4, 1 .4). Error bars denote 
standard deviations in wavenumber bins. The best-fit lines and slopes, which 
were obtained over the range 3 < k < 35, are also shown. 
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Figure 2. Spectral slopes (7) of lognormal density fields vs. the original 
spectral slopes (7 g ) of the Gaussian density fields. Best-fit polynomial curves 
for every density variance cr\ np are also shown. The circles and curves were 
shifted vertically by the amounts of A for clarity. 



density variance-Mach number relationship, using MHD sim- 
ulations. The relationship closely follows the form of the 
3D density variance-Mach number relationship but includes 
a scaling parameter A such that 

ai N =Aln(l+b 2 M 2 J, (1) 

where a 2 nN is the variance of the logarithm of the column den- 
sity N 4 . 

In the present study, quantitative expressions relating the 
3D density variance, spectral slope, and column density vari- 
ance are obtained by using the fractional Brownian motion 
(ffim) algorithm (Saupe 1988; Stutzki et al. 1998; Elmegreen 
2002). A lognormal density field is obtained by generating a 
Gaussian random field with a power-law spectral slope (j g ) 
and exponentiating the field (Elmegreen 2002). This letter is 
organized as follows. In Section 2, we describe how to infer 
the variance of the 3D density field from the given spectral 
slope and column density variance. Section 3 investigates the 
Ofaff— M s relationship. Concluding remarks are given in Sec- 
tion 4. 

2. MODEL OF THE LOGNORMAL DENSITY FIELD 

The fBm structures are generated by first assigning 3D 
Fourier coefficients following the prescription of Elmegreen 
(2002). The Fourier amplitudes are generated to be distributed 

as a normal Gaussian with the variance < |a(k)| 2 >= |k| 7s . 
The inverse Fourier transform gives a Gaussian random field 
p g (x). We then multiply the density field with the desired stan- 
dard deviation of the logarithmic density <r\ np , and exponenti- 
ate the field to obtain p = exp(a-\ np p g ). 

The standard deviation <t\„n depends on both 7 g and <7i np . 
We should also note that the spectral slopes 7 of resulting log- 
normal density fields are different from the input slopes 7 g . 
Therefore, we need to find the relationship between the input 
variables (7 g , <j\ Dp ) and the outputs (7, <Ti n jv). In addition, the 
generated density fields with the same parameters (7 g , <j\ np ) 

4 Note that o- lnN / N(j = (t\ u n, but cr N / No = on /No for the mean column den- 
sity No ■ The same relationships are applied to the 3D density p. 
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Figure 3. The column density deviations o\ nN and spectral slopes 7 of log- 
normal density fields obtained for various combinations of tri n „ and 7 g of the 
Gaussian random field. Solid and dotted lines are best-fit curves tracing the 
(cri n jv, 7) pairs corresponding to constant oi Rp and jg, respectively. 

show large fluctuations due to random phases and amplitudes 
in Fourier space. We thus generated a large volume of log- 
normal density fields by varying 7 g and a\ n p and obtained the 
average relationships between the parameters. We varied 7 g 
from 2.0 to 5.25 in steps of 0.25 and <j\ np from 0.01 to 2.76 
in steps of 0.25. One hundred random realizations with a box 
size of 128 3 were generated for each combination of 7 g and 
cr\a P , resulting in a total of 16,800 realizations. We also gen- 
erated 30 realizations with a resolution of 256 3 for every (7 g , 
cinp) pair and confirmed that the 256 3 dataset is consistent 
with the results obtained with the 128 3 resolution. As will 
be shown later, the present results are also consistent with 
the limited samples obtained with a box size of 512 3 . The 
standard deviations and spectral slopes were calculated for 
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Table 1 

Polynomial coefficients for 7 as a function of (7g, o"i np ), and for 7 g and 
<ti„£ as functions of (7, <Ti np ). 
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Figure 4. Column density variance-Mach number relationship. The up- 
per three black-curves represent the cases of compressive (b = 1), mixing 
(b = 0.5), and solenoidal mode (b = 1/3). The lowest black curves are ob- 
tained by assuming b = 0.29 and 7 = 3.81M.7 16 -0.3 to mimic the magnetic 
field effects. Dashed black lines represent the best-fit results to Equation ( 1 ) 
with the parameters (b, A) denoted along the curves. Dashed blue line shows 
the curve of Burkhart & Lazarian (2012). Forcing types (c, m, s) and the sim- 
ulation regimes (H, M) are also denoted. Red lines are obtained by varying 
the 7 -M s relationship arbitrarily and demonstrate that Equation (1) is not 
satisfied unless proper relationships between 7 and M s are assumed. 

density fields with a shallower spectral slope. For instance, 
the density fields with 7 = have no spatial correlation. 

The equations relating 7 and 7 g will be useful in generat- 
ing the realistic ISM density structures. The equation relating 
cinAr to 7 and <j\ np or Figure 3 provides a tool to derive <j\ np 
from the observable a\ n n and 7. Once we obtain 7 and a\ a ^ 
from observations, we can vary a\ np to find the best-fit value 
that gives the observed a\„N- 

3. COLUMN DENSITY VARIANCE-MACH NUMBER RELATIONSHIP 

We now investigate the relationship between cri n Ai and M s 
by using the expressions derived in Section 2. In the above, 
we have assumed that 7 and <j\ np are independent variables. 
However, 7 and cr\ np are not independent in the real turbulent 
media, but rather depend on M s . If we ignore magnetic fields, 
the statistical properties are uniquely determined by M s and b. 
We therefore need to specify the dependencies of 7 and a\ n p 
onM s . 

Seon (2009) combined the results of Kim & Ryu (2005), 
Kritsuk et al. (2006), and Padoan et al. (2004), and obtained a 
simple relation between M s and 7, i.e., 7 = 3.81M" 016 , which 
may be applicable to solenoidal forcing. The power spectra of 
the compressively driven turbulence are considerably steeper 
than those of the solenoidal forcing. The spectral slopes at 
M s pa 2.3 and 5.6 for compressive forcing are provided by 
Schmidt et al. (2009) and Federrath et al. (2009), respectively, 
and are ~ 0.6 larger than the values for the solenoidal case. 
We therefore assume that 7 = 3.81M; 16 + 0.6 for the com- 
pressive mode. There are two transverse (solenoidal) and 
one longitudinal (compressive) wave mode in the 3D space. 
Therefore, the solenoidal and compressive forcing would be 
mixed with a ratio of 2: 1 if no preferential forcing mode were 
provided (Federrath et al. 2010). The power spectrum of this 
natural mixing mode might be obtained by adding two power 
spectra for solenoidal and compressive forcing with a ratio 
of 2: 1 . The resulting power spectrum is also represented by 
a power-law with a spectral slope, which is at most ~ 0.06 
larger than the slope of the solenoidal case. We therefore as- 



Note: The numbers in parentheses are exponents on 10, meaning that 
6.609(-l) = 6.609x10-'. 

every realization and averaged to obtain the ensemble aver- 
age values of (7, o"i n w) for each combination of the input (7 g , 
<j] np ) values. The power spectra of resulting lognormal den- 
sity fields show a slight flattening at large wavenumbers. Ex- 
amples of the spectra are shown in Figure 1 . Spectral slopes 
were obtained by least-squares fits over a wavenumber range 
of 3 < k < 35. Here, the dimensionless wavenumber is de- 
fined by k = L/X with the wavelength A and box size L. 

For ease of use, we parameterized the relationships between 
the parameters 7, 7 g , <J\ np , and a\ n ^ using polynomials. For 
every a\ ap , a\ a N was parameterized with a cubic function of 
7, i.e., 

3 

C]aN = 22^' =C + C 1 7 + C27 2 + C37 3 . (2) 

;=o 

The coefficients c, (i = 0, . . . , 3) were then fitted with a quartic 
polynomial of a\ np , i.e., 

4 

CO = X^Wnp = P0 + Pl&\np + P2<?l p + P3<?l np + P4<?ip- O) 

In this way, a\ n N could be estimated as a function of 7 g and 
(Ti n p . We also parameterized 7 with a cubic function of 7 g 
and 7 g with another cubic function of 7 for every <j\ np , i.e., 
7 = S,a,7 g and 7 g = £,'fr,7 ! (i = 0, ... ,3). The coefficients a,- 
and bj (/ = 0, . . . , 3) were then fitted with quartic polynomial 
functions of <j\ np , as done for the coefficients c,-. 

The final coefficients are shown in Table 1 . Figure 2 shows 
the resulting 7's as functions of 7 g for every <r\ np . Figure 3 
presents the resulting (7, crin/v) pairs for every (7 g , <j\ np ) pair. 
In the figures, the solid and dotted lines are the best-fit poly- 
nomials to reproduce the (7, erin/vO pairs corresponding to con- 
stant (j\ np and 7 g , respectively. 

The power spectra of lognormal density {p) fields are usu- 
ally shallower than those of the original Gaussian density 
(Inp) fields, i.e., 7 < j g . The property that the power spec- 
tra of Inp are stiffer than those of p in turbulence simula- 
tions is also noticeable in Table 2 of Kowal et al. (2007). The 
difference between 7 and 7 g becomes more significant for a 
larger <r\ np and for shallower (smaller) 7 g . Column density 
variance is always smaller than the variance of 3D density, 
as expected. The ratio of a^N to a\ np becomes smaller for a 
shallower 7. This is caused by less spatial correlation in the 
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sume the same spectral slope as the solenoidal forcing case 
for the natural mixing case. 

The density variance in logarithmic scale is obtained from 
the variance in linear scale using the property of a lognor- 
mal distribution, i.e., af np = \n{l + a 2 p , ) = ln(l + b 2 M;). The 
forcing parameter b is 1/3, 0.5, and 1 for solenoidal, natural 
mixing, and compressive modes, respectively (Federrath et al. 
2008, 2009, 2010). 

Using 7 and cr\ np as functions of M s , and Equation (2), we 
calculate a\ n n for solenoidal, natural mixing, and compressive 
modes as functions of M s . The results are shown in Figure 4, 
wherein the upper three black-curves show <j\ n N as functions 
of M s for three different forcing types. We fitted the obtained 
curves with Equation (1). The best-fit A values are found to 
be 0.2, 0.24, and 0.38 for the solenoidal, natural mixing, and 
compressive modes, respectively. In Figure 4, the b and best- 
fit A values are denoted in parentheses along the curves. The 
best-fit curves are shown as dashed lines. 

Applying the expression relating 7 and 7„ derived with the 
128 3 dataset, we also generated the lognormal density fields 
with a higher resolution of 512 3 for M s = 0.5, 1, 2.5, 5, 10, 
and 20. For every combination of parameters (b, M s ), we pro- 
duced 30 realizations and calculated the ensemble average of 
tJiniv and its statistical dispersion. Circles and error bars in 
Figure 4 represent the ensemble averages and the dispersions, 
respectively, which closely follow the curves calculated with 
the 128 3 dataset. We therefore conclude that the expressions 
obtained with a rather lower resolution of 128 3 are also ap- 
plicable to the higher resolution of 512 3 . The error bars ob- 
tained with the 5 12 3 resolution indicate the ranges of variation 
between realizations. 

In Figure 4, we over-plotted the data points obtained 
from the HD simulations with solenoidal (blue squares) and 
compressive forcing (red squares) modes of Federrath et al. 
(2010), who provide the standard deviations of linear column 
density (ctjv/jv ) and logarithmic column density (cin^). The 
filled squares show the <j\ nN values that they directly measured 
with the logarithmic column density and the open squares rep- 
resent the values [ln(l + <T ^/ Wo )] 2 calculated with the stan- 
dard deviations of the column density a N / Na assuming that the 
column density PDF is a perfect lognormal distribution. The 
values estimated with the two different methods show a sig- 
nificant discrepancy, especially for compressive mode, indi- 
cating that the column density PDF is not a perfect lognormal 
(see also Price et al. 2011 for the discrepancies). Indeed, the 
PDFs for compressive forcing show large departures from a 
perfect lognormal function (Federrath et al. 2010). Even with 
this departure from the perfect lognormal function, the present 
results are consistent with their results within the variation 
ranges denoted by error bars in the figure. 

In the above analyses, we have assumed HD cases. 
Burkhart & Lazarian (2012) fitted the MHD simulation re- 
sults for solenoidal forcing with Equation (1) and obtained 
a best-fit A of 0.11, which is different from our result of 
A = 0.2. The difference might be attributable to the mag- 
netic field effects. First, the density variance in magnetized 
gas is significantly lower than that in the HD counterparts 
(Ostriker et al. 2001; Price et al. 2011; Molina et al. 2012). 
The effect could be mimicked by reducing b below the value 
of the HD case, while assuming the same relationship between 
7 and M s for the HD case. We found that the of nA ,-M s re- 
lationship derived with b = 0.22 is very close to the result of 



Burkhart & Lazarian (2012). Second, the spectral slope of the 
strongly magnetized media appears to be shallower than those 
of unmagnetized or weakly magnetized media. Padoan et al. 
(2004) found that the MHD model with approximate equipar- 
tition of kinetic and magnetic energies yields a shallower den- 
sity power spectrum than the super-Alfvenic model. The re- 
sults of Kowal et al. (2007) also provide some indications that 
the power-law slopes of the supersonic, sub-Alfvenic mod- 
els are slightly smaller than those of super-Alfvenic models. 
Their spectral slopes for the cases of a strong magnetic field 
are lower than the values expected from 7 = 3.81M" 016 by 
amounts of ~ 0.3-0.4. Indeed, we could reproduce the best- 
fit curve of Burkhart & Lazarian (2012) with a spectral slope 
of 7 = 3.81M-° 16 -0.5 and b = 1/3. In a real situation, both 
effects may play roles together. Therefore, we lowered both 
7 and b and found that the results with 7 = 3.81M~ 016 -0.3 
and b = 0.29 match the result of Burkhart & Lazarian (2012), 
as shown by the lowest black curve in Figure 4. In the figure, 
the black and blue dashed lines denote the best-fit results to 
Equation (1) with b = 0.29 and b = 1/3, respectively. 

Burkhart et al. (2010) presented 7 as a function of M s for 
MHD simulations with different magnetic field strengths and 
found that the magnetic field effect on 7 is more important at 
lower M s than at higher M s . Therefore, we may need further 
investigations on the magnetic field effects. 

Another aspect worth noting is that we can constrain the 
condition for Equation (1) to hold, by varying the equation 
between 7 and M s . For instance, in Figure 4, we plotted the 
eri n Ai-M s relationship by assuming a constant spectral slope 
of 7 = 2.8 (upper red curve) for b = 0.5. It is obvious that the 
resulting curve cannot be expressed by Equation (1). As we 
assume 7 increases with M s , we obtain more rapidly increas- 
ing curves than the red curve. On the other hand, if we assume 
7 decreases more rapidly than 7 = 3.81M" 016 with increasing 
M s , the resulting curve shows a much slower increase with 
increasing M s , as indicated by the lower red curve obtained 
for 7 oc M" 35 . Therefore, we conclude that the necessary 
condition for Equation (1) is that the spectral slope should 
moderately decrease as M s increases. 

4. CONCLUDING REMARKS 

We investigated the dependency of the column density vari- 
ance on the 3D density variance and the power-law slope 
of the density power spectrum by using the fBm structures. 
Adopting the spectral slopes and 3D density variances that 
are appropriate in the HD regime, we obtained the same rela- 
tionship between the column density variance and the Mach 
number as suggested by Burkhart & Lazarian (2012), but with 
a different scaling parameter A. However, when the magnetic 
field effects were included, the parameter A agreed with the 
MHD result of Burkhart & Lazarian (2012). Magnetic fields 
prevent the turbulent gas from reaching very low densities as 
well as very high densities (Molina et al. 2012), and conse- 
quently have important implications for the models of the star 
formation rate (e.g., Federrath & Klessen 2012). 

The approach adopted in this letter is similar to that 
of Padoan et al. (1997) in the point that they also use 
random density fields to derive the 3D density variance. 
The fBm structures also have been used to investigate the 
density structure of molecular clouds (Stutzki et al. 1998) 
and the stellar initial mass function (Elmegreen 2002; 
Shadmehri & Elmegreen 2012). Note that Stutzki et al. 
(1998) neglected the dispersion in Fourier amplitudes in gen- 
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erating Gaussian random fields. We compared the cases of 
fixed Fourier amplitudes with the results presented here, and 
found that the fixed amplitudes give rise to slightly stiffer 
spectral slopes than the present results. 

The present results can be applied directly to the obser- 
vations to infer the 3D density variance from the observed 
column density variance and spectral slope. We can vary 
o"inp in Equations (2) and (3) to find the best-fit value that 
matches the observed <j\ d n- In applying the results, we need 
to consider the effects of instrument noise and telescope 
smoothing. However, the inclusion of noise and smoothing 
may not severely affect the applicability of the results, as in 
Burkhart & Lazarian (2012). 

The density fluctuation of the ISM plays a crucial role in un- 
derstanding the propagation of radiation. The intensity PDF 
of the dust-scattered starlight in far-ultraviolet was found to 
be lognormal, which might be caused by the density struc- 
ture of interstellar dust (Seon et al. 201 1). Escape of the ion- 
izing radiation field through the clumpy ISM has been ar- 
gued to be the origin of ubiquitous diffuse ionized gas (e.g., 
Haffner et al. 2009). This has not, however, been clearly con- 
firmed using realistic ISM density structures (c.f., Seon 2009; 
Seon & Witt 2012). The present study provides a practical 
method to generate the realistic density structures for the ra- 
diative transfer problems. 
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